Magnetostatics solution (linear case)

%%capture
%run current_density.ipynb

##############################################################################
# Tree/Cotree gauging
##############################################################################
MESH = pde.mesh3.netgen(geoOCCmesh)
R = pde.tools.tree_cotree_gauge(MESH)

R

<27728x23597 sparse matrix of type '<class 'numpy.float64'>'
        with 23597 stored elements in Compressed Sparse Column format>

##############################################################################
# Assembly
##############################################################################

order = 1

phi_H1 = pde.h1.assemble3(MESH, space = 'P1', matrix = 'M', order = order)
dphix_H1, dphiy_H1, dphiz_H1 = pde.h1.assemble3(MESH, space = 'P1', matrix = 'K', order = order)
D = pde.int.assemble3(MESH, order = order)
# R0, RSS = pde.h1.assembleR3(MESH, space = 'P1', faces = 'ambient_face')

M = phi_H1 @ D @ phi_H1.T

K = dphix_H1 @ D @ dphix_H1.T +\
    dphiy_H1 @ D @ dphiy_H1.T +\
    dphiz_H1 @ D @ dphiz_H1.T

# Kn = RSS @ K @ RSS.T

phix_Hcurl, phiy_Hcurl, phiz_Hcurl = pde.hcurl.assemble3(MESH, space = 'N0', matrix = 'M', order = order)
curlphix_Hcurl, curlphiy_Hcurl, curlphiz_Hcurl = pde.hcurl.assemble3(MESH, space = 'N0', matrix = 'K', order = order)

M_Hcurl = phix_Hcurl @ D @ phix_Hcurl.T +\
          phiy_Hcurl @ D @ phiy_Hcurl.T +\
          phiz_Hcurl @ D @ phiz_Hcurl.T

K_Hcurl = curlphix_Hcurl @ D @ curlphix_Hcurl.T +\
          curlphiy_Hcurl @ D @ curlphiy_Hcurl.T +\
          curlphiz_Hcurl @ D @ curlphiz_Hcurl.T

C_Hcurl_H1 = phix_Hcurl @ D @ dphix_H1.T +\
             phiy_Hcurl @ D @ dphiy_H1.T +\
             phiz_Hcurl @ D @ dphiz_H1.T

curlphix_Hcurl_P0, curlphiy_Hcurl_P0, curlphiz_Hcurl_P0 = pde.hcurl.assemble3(MESH, space = 'N0', matrix = 'K', order = 0)
phix_Hcurl_P0, phiy_Hcurl_P0, phiz_Hcurl_P0 = pde.hcurl.assemble3(MESH, space = 'N0', matrix = 'M', order = 0)
dphix_H1_P0, dphiy_H1_P0, dphiz_H1_P0 = pde.h1.assemble3(MESH, space = 'P1', matrix = 'K', order = 0)

KR = R.T@K_Hcurl@R
MR = R.T@M_Hcurl@R

r = jx_hdiv @ D @ phix_Hcurl.T +\
    jy_hdiv @ D @ phiy_Hcurl.T +\
    jz_hdiv @ D @ phiz_Hcurl.T

# # use this instead if you use the H1-L2 approx for the current J
# r = jx @ D @ phix_Hcurl.T +\
#     jy @ D @ phiy_Hcurl.T +\
#     jz @ D @ phiz_Hcurl.T

cholKR = chol(KR)
a = cholKR.solve_A(R.T@r)
a = R@a

Bx = curlphix_Hcurl_P0.T @ a
By = curlphiy_Hcurl_P0.T @ a
Bz = curlphiz_Hcurl_P0.T @ a

C:\Users\Radu\AppData\Local\Temp\ipykernel_15748\2036867034.py:51: CholmodTypeConversionWarning:

converting matrix of class csr_matrix to CSC format

##############################################################################
# Storing to vtk
##############################################################################

grid = pde.tools.vtklib.createVTK(MESH)
pde.tools.vtklib.add_L2_Vector(grid,Bx,By,Bz,'grad_x')
pde.tools.vtklib.writeVTK(grid, 'magnetostatics_solution.vtu')

import pyvista as pv
mesh = pv.read('magnetostatics_solution.vtu')
mesh2 = pv.read('current_density.vtu')

mesh.set_active_scalars("Scalars_")
threshed = mesh.threshold([0,1])

p = pv.Plotter()
p.add_mesh(threshed, style='surface', color = "w", opacity=0.2, label=None)

threshed.set_active_vectors("grad_x")
arrows = mesh.glyph(scale="grad_x", orient=True, tolerance=0.03, factor=1000)
p.add_mesh(arrows, color="orange")


mesh2.set_active_vectors("J_L2")
arrows2 = mesh2.glyph(scale="J_L2", orient=True, tolerance=0.03, factor=9500.0)
p.add_mesh(arrows2, color="black")

p.camera_position = [(0, -500, 200),(0, 0, 0),(0, 0, 0)]
p.show(jupyter_backend="html")